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ABSTRACT 


- Investigation of steady and unsteady flowfields over airfoils is an active 
area of current computational and experimental research. In this study, the 
compressible, viscous flow over a single and multi-element airfoil is numerically 
simulated by solving the Navier Stokes equations. The motivation for this work 
includes interest in studying the effects of a stationary/flapping airfoil combination 
in tandem configuration. A single-block Navier- Stokes (NS) solver is employed 
to compute unsteady flowfields. Turbulence is treated using the Baldwin-Lomax 
turbulence model. A single C-grid is generated and it is partially distorted to 
simulate the flapping motion. Numerical solutions are obtained for flows at a fixed 
angle of attack and for unsteady flows over flapping airfoils. The numerical 
solutions agree well with the experimental data. The difficulties faced during the 


study are discussed and future improvements are suggested. 
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L INTRODUCTION 


Studies of steady and unsteady flowfields over multi-element airfoils is an active 
area of numerical and experimental research. While notable progress is being made in the 
development of computing methods to analyze flows around single airfoils, methods to 
compute more complex geometric configurations are now being developed. In this study, 
a Navier Stokes solver was modified to accommodate the complexity of a multi-element 
configuration to study the thrust generating effect of a stationary /flapping airfoil 
combination. 

Thrust generation due to airfoil flapping was recognized by early researchers, such 
as Knoller [Ref. 1], in explaining the bird’s ability to generate a propulsive force by means 
of flapping their wing. Previous experimental investigations of unsteady flows over 


oscillating airfoils by Neace [Ref. 2] ,Tuncer and Platzer [Ref. 3], and Dohring [Ref. 4] 





have shown a propulsive effect and the associated efficiencies. Figure 1 shows how the 
airfoil motion creates an induced velocity component. This allows for the generation of a 


forward or thrust component of force. 


Thrust 





Figure 1 Propulsive Force on Plunging Airfoil. Neace [Ref. 2 ] 


In the process of this research, a single block Navier Stokes solver was first 
validated on a Ames-01 airfoil. The solver was subsequently modified to handle multi- 


element airfoil configurations utilizing a single structured C-grid. Computational grids 





were generated utilizing GRIDGEN . Other grid generating tools were evaluated, but 
only GRIDGEN provided the flexibility to properly distribute the gridlines in critical areas. 
The initial grid configuration utilizing NACA 4412 and 4415 airfoils within a wind tunnel 
boundary was used in the validation process of the code. Two other configurations were 
studied under unsteady conditions . In the unsteady flow case the leading airfoil is kept 
stationary while the trailing airfoil undergoes a flapping motion. In Case I, flap is placed 
at .Sc aft of the main airfoil and in Case II at .25c aft of the main airfoil These two cases 


were studied at different flapping frequencies and amplitudes. 





i. METHOD 
A. GRID GENERATION 
1; Code Validation Grid 


During the grid generation process proper distribution of grid points around the 
airfoil’s leading edge and trailing edge, orthogonality of the grid lines on the surface, and 
smooth variations of the grid density are the most critical factors. In order to accurately 
calculate the flow gradients in the direction normal to the surface in the boundary layer, it 
is necessary to make the normal grid spacing very fine close to the solid surfaces. Single 
C-type grids for NACA 4412 (Main airfoil) and NACA 4415 (Flap) were generated for 
the code validation process modeling the configuration in Figure 2. The grids included in 
this report were generated using GRIDGEN software packages. 

Adair and Horne [Ref. 5] conducted experiments with this configuration in the 7 
by. 10 Foot Wind Tunnel at NASA Ames Research Center. In the experiment, the chord 
length of the main airfoil is .90 m and that of the flap is 36m. The geometric location of 
the flap relative to the main airfoil was described by the flap gap (FG) , the flap overlap 
(FO) and the flap deflection (5,). For this case FG = 0.035c¢ , FO = 0.028c , d¢ = 21.8° 


and the main airfoil angle of attack («) was set to 8.2°. 
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Figure 2 Experimental Configuration Adair and Horne [Ref.5] 
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The proximity of the airfoils, and the sensitivity of the solver puts a heavy 
constraint on the computational grid around the downstream airfoil. In this case it was 
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computationally demanding and the most stable. This method produces the smoothest 
Manual [Ref. 6]. See sample of GRIDGEN input file and elliptical solver parameters in 


orthogonality. Among those, the Laplace control function with a .7 relaxation parameter 
possible grid. For a detailed description of all related parameters see the GRIDGEN 


was the standard selection for this research. This type of function 1s the least 


the Appendix B. 








The following is a list of parameters and procedures used to run the elliptical solver: 
e Relaxation factor = 0.7 
e Laplace control function 
e Set BC’s on the airfoil surface to orthogonality. By default all points on the 
edges will remain fixed as the elliptical solver is run. 

e Turn foreground BC’s on (elliptical solver menu) with default values. 

e Run the solver for critical areas separately, this will expedite the process and 

will result in a smoother grid. 

e Run the solver until an acceptable grid is achieved. 

Figure 4 shows the grid distribution in the wake area for the initial configuration 
and the up-flap configuration. For the aft-flap configuration the flap was moved aft 0.03c, 
for the up-flap configuration the flap was moved up from the aft position 0.03c. Due to 
the positioning of the flap in the wake of the main airfoil, (up-flap configuration) the 
density of the grid lines between airfoils was increased. The grid dimensions for the initial 


and up-flap configurations were 361 x 100 and 463 x 100, respectively. 
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Figure 5 shows the full wind tunnel model for the aft-flap configuration used 


during the validation process. 
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Figure 5 Aft-flap Wind Tunnel Grid 


2. Unsteady Case I and II Grids 


For unsteady computations two separate grids were generated as Case I and II. 
NACA 4412 and 4415 airfoils were used for these computations. The airfoil combination 
was set in a tandem configuration with a gap between airfoils of .5c for Case I and .25c 
for Case II. (Figure 6). The length of the flap (NACA 4415) was set to 0.33c. The main 
airfoil was set at « = 10° and the flap at zero angle of attack with respect to the main 
airfoil chord line. The grid dimensions are 391 x 100 for Case I and 371 x 100 for Case 
IL. 

These configurations eliminate some of the complexities of the initial configuration 
grid allowing the evaluation of the oscillation effects without the influence of any other 
variables. In addition, having the trailing airfoil far enough from the main airfoil facilitates 
the implementation of distortion required to model the oscillation process. Details of the 


erid distortion are discussed in the next section. 
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2-D, thin layer Navier- 


Figure 6 Case I (top) and II grids(bottom) 


thin layer, Navier-Stokes solver with the third order accurate Osher's 
is given as follows 


> 


Navier-Stokes Solver 


ALGORITHM 
An implicit 


L. 


upwind biased flux difference splitting scheme is employed to compute the flowfield in the 
Stokes equations in a curvilinear coordinate system, (€,¢), along the axial and normal 


computational grid. The strong conservation law form of the 


direction respectively, 


B. 
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S is the thin layer approximation of the viscous fluxes in the ¢ 


Q is the vector of conservative variables 


inviscid flux vectors, and 


Where 


direction normal to the airfoil surface: 








pu pw 
iA) POT te? Af) Pee 
S| pwutEp | > S| pwW +p 
(e+ p)lU-¢.p (e+ pW-¢,p 
0 


1 pmyu, +(e /3)m,¢, (2.2) 
J MMW, +(u/3)m,¢, 

pm, + (u6/3)m, +(G,utow 
where | 
m=Co+0 (2.3) 
Mm, =6,U,+¢,w, 


4{ & 
m, =(u° +w*)/2+KPr ies 


and U ,W are the contravariant velocity components. In the above equations, all 
dimensions are normalized with the airfoil chord length, c. p is the density normalized 
with the free-stream density ,; u and w are the Cartesian velocity components in the 


physical domain, which are normalized with the free-stream speed of sound a. e is the 


total energy per unit volume normalized with p,4.,: and Pr is the Prandtl number. The 


pressure is related to density and total energy through the equation of state for an ideal 


gas, p= (vy —le— plu? +w’)/2]. (2.4) 

The flowfield is assumed to be fully turbulent. The turbulence modeling is used to relate 
the Reynolds shear stress to the local mean-velocity gradient allowing numerical flow field 
calculation. The Baldwin--Lomax algebraic turbulence model is currently implemented. 
This model is a two layer eddy viscosity model which simulates the effect of turbulence in 
terms of eddy viscosity coefficient. A complete description of the model is given in 


Baldwin [Ref. 7]. 














2. Boundary Conditions 


The computational domain includes both airfoil surfaces and extends ten chord 
lengths away from the airfoils. Boundary conditions are applied on the airfoil surfaces, and 
at the farfield boundaries. On the airfoil surfaces the no slip boundary condition is 
applied. For the flapping airfoil, the surface fluid velocity is set equal to the prescribed 
airfoil flapping velocity so that the no slip condition is satisfied. Since the formulation of 
the Navier-Stokes solver is based on an inertial frame of reference, the flapping motion of 
the airfoil is implemented by moving the airfoil and the computational grid around it in the 
transverse direction as described by the frequency and the amplitude of the flapping 
motion by: | 


h = —Acos(at) (2.5) 


Where A denotes the mean amplitude normalized with the chord length. is the 
frequency of the oscillatory flapping motion, which is used in terms of reduced frequency, 
k=ac/2U, (2.6) 
At the farfield inflow and outflow boundaries the flow variables are evaluated 
using the zero order Riemann invariant extrapolation. 
For the code validation studies, no slip boundary condition is applied at the tunnel 
walls. 


3, Unsteady-Motion Routine 


The features pertinent to the grid distortion in the NS code were modified so that 
the trailing airfoil can move in the cross-flow direction with respect to the leading airfoil. 
During the flapping motion the main airfoil remains stationary while the flap oscillates in 
the cross flow direction. This motion is introduced by a partial grid distortion. The wake 
region between airfoils is the area of concern since the distortion is gradually introduced 
from the trailing edge of the main airfoil to the flap leading edge. During this distortion 
the integrity of the grid in terms of orthogonality, smoothness and clustering of grid point 
should remain physically valid for a solution to converge. In order to maintain a smooth 


variation of grid distortion a hyperbolic tangent is used as the transition function between 
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the stationary and moving part of the grid.. The grid’s distortion adjusts smoothly using 
this function. Figure 7 shows this function. The grid motion routine (GRMOVE 
subroutine of the NS code) is given in the Appendix A. 


TE Main 
Airfoil 


Position at maximum 
distortion. h=-A 





Wake Area 


Figure 7 Transition Function- Hyperbolic Tangent 
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Iii. RESULTS AND DISCUSSION 


The validity of the code was first tested with a specified configuration. The code 
was tested under the steady-state conditions described below. The solution was then 
evaluated for a number of configurations and compared to the experimental data. 
Following the validation process the flowfield with the flapping trailing airfoil was studied 
for two different configurations. 


A. CODE VALIDATION 


1. Conditions--Steady State 


For the validation process the flowfield was computed at: 

°ea=82° 

e Mach= .2 

e Reynolds number 1.8 x 10°, fully turbulent. 
The solver was first run with the initial wing-flap configuration. Convergence was 
reached at 1500 iterations based on the variations in the flow residuals and the 
aerodynamic loads. Similarly, the solver was run for the aft-flap and up-aft flap 
configurations. A typical NS input file is shown in the Appendix A. Figure 8 shows the 


convergence history of the initial case and Figure 9 shows the aerodynamic loads for the 


same configuration. 


Residual 





1000 1500 2500 


iterations 


2000 3000 


Figure 8 Convergence History Initial Configuration 
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Figure 9 Aerodynamic Loads - Initial Configuration 


Figure 10 shows the pressure coefficient distribution for all cases. The pressure changes 


indicated in this figure for the different cases shows a general trend that represents a 


e@ Exp Data 
x Aft-Flap 
Oo Up-Flap 
| alnitial Configuration 





Figure 10 Computed vs. Experimental Cp Distribution 
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good agreement with the experimental data. Figure 11 summarizes the effects of the gap 
on the separated flap flow. This figure shows how the flowfield is affected as the flap is 
lowered or moved aft. As the gap is increased more fluid from below the main airfoil 
moves into the flap through the gap pushing the flow farther away from the flap leading 
edge, Alemdaroglu [Ref. 8]. As the gap is increased the flow will move the streamlines 
upward affecting the flow separation over the flap. This condition has significant effect on 
the pressure distribution over the aft-upper surface of the main airfoil as noted in F igure 
10. 


Mass -Flux 








Figure 11 Computed flowfield - Initial Configuration (top), Aft-flap at b= 13.6° 
(bottom) 
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B. UNSTEADY CASE 


This part of the research focused on the aerodynamic interaction between two 
NACA airfoils as shown in Figure 6. The main airfoil (NACA 4412) is stationary while 
the flap (NACA 4415) undergoes a flapping motion. The flowfields were evaluated under 
the following conditions: 

e «= 10 degrees 

e Mach = 3 

e Reynolds number 3.93 x 10° 
The unsteady computations were initialized with a steady-state solution computed at the 
maximum flapping amplitude where the flapping velocity is zero, Tuncer and Platzer [Ref. 
| 3]. See a copy of the input file in the Appendix A. The unsteady computations were 
carried out for up to three periods of harmonic flapping motion. The flowfields presented 
were taken from the last period of the computation cycle. All unsteady runs were 
completed at the Naval Postgraduate School Cray Y-MP J94. The average computer 


processing time (CPU) for 15,000 time step of unsteady computations was 9.21 hours. 


1. Case I 


The flowfield for this case was initially computed with a flapping motion of A= 
.10c and k= 0.5. The computed streamlines for both the steady and unsteady case are 
shown in Figure 12. It can be noted in this figure that the fluid particles no longer travel 
freely through the gap to the upper flap section. Instead, the suction created by the 
flapping motion keeps the dividing streamline closer to the flap leading edge. The flow 


fields for both the steady-state and unsteady solution are shown in Figure 13. 
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Figure 13 Velocity Magnitude Computed for Steady-State (top) and Unsteady (bottom), 
A=0.10c, k=.5 


Due to the gap size the effect on pressure distribution over the main airfoil was very small 
as shown in Figure 14. On the other hand, the boundary layer profile changes, shown in 
Figure 15, at the aft upper surface of the main airfoil are significant. It is noted that the 
flow reversal is reduced , the magnitude of the velocity vectors is increased and the 
boundary layer reattaches to the surface close to the trailing edge. The changes are really 
noticeable in the trailing edge region which indicates that the suction produced by the 
flapping airfoil is not sufficient for a reattachment of the boundary layer forward of that 


region. 
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Figure 15 Boundary Layer Profile, Steady State (top), Unsteady (bottom), 
A =0.10c, k =.5 














The case I configuration was tested at A = 0.05 and k = .5, with very similar 
results, as shown in Figure 16. This figure suggests that due to the reduced oscillating 


amplitude the suction effect was significantly reduced. 





Figure 16 Velocity Magnitude Computed for Steady-State (top) and Unsteady (bottom), 
A=0.05c, k=.5 


The time history of the unsteady aerodynamic loads of the flapping motion is 
shown in Figure 17. Following the initial transition , all aerodynamic loads attain a 
periodic behavior with respect to the flapping motion. For the lift coefficient the 
oscillation created peak variations of less than one percent for both cases. Similarly, drag 


coefficient variations were less than 2.5 percent for A=.1c and one percent for A=.05c. 
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Figure 17 Time History of the Unsteady Lift and Drag Coefficients 


Figures 18 and 19 shows the particle traces for four different flap positions. 
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Figure 18 Particle Traces, Lower Position (top), Center Position (bottom) 
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Figure 19 Particle Traces , Traveling Up through Wake (top), Upper Position (bottom) 
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2. Case Il 


For this case the gap between airfoils was changed to .25c. Solutions were 
obtained for A= .05 and k =.5 and 1.0. Figure 20 shows how changes in oscillating 
frequency influence the thrust generation. As frequency increases, the suction increases 


pulling the wake of the main airfoil closer to the flap boundary. 
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Figure 20 Velocity Magnitude Computed for Steady-State (top), k = .5 (center) and 
k=1.0 (bottom) 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


Subsonic flow over stationary airfoil/oscillating flap combinations were computed 
using a Navier-Stokes solver. The computed flowfields confirm the observations made in 
recent experiments performed by Dohring [Ref. 4]. However additional computations for 
other configurations at different frequencies, amplitudes, angle-of-attack, airfoil sections, 
etc., are required before more quantitative conclusion can be drawn. Preliminary review 
of the computed boundary layer characteristics indicates that with a proper configuration 
and applying a proper combination of the flapping parameters mentioned above a 


reattachment of the boundary layer may result. 


B. RECOMMENDATIONS 


To fully satisfy the goals of this research a wider parametric study is required. In 
order to obtain a clear picture of the flapping effects and to be able to quantify the thrust 
efficiency the effect of the following parameters listed below needs to be identified. 

e Frequency of oscillation 

e Amplitude | 

e Reynolds number 

e Gap size and flap location 

e Angle-of-attack 

e Airfoil section (i.e. symmetric, non-symmetric, flat plate, etc. ) 

In addition; | 

1. The complexity of the flow, especially in the wake area of both airfoils, 
requires special attention. For complex turbulent flows, the choice of turbulence model 
will have a significant impact on accuracy and computational time. Therefore, other 
turbulence model, such as k-e or Baldwin-Barth, needs to be investigated for the multi- 


element configuration, Nelson [Ref. 9]. 
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2. Grid generation is critical toward obtaining accurate solutions. Tools currently 
available are complex. Therefore, more work is required to create a grid that solves the 
overhang and overlap problem discussed in the grid generation section. A multi-block 


grid may solve this problem. 
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APPENDIX A 


A. DESCRIPTION OF NS INPUT FILE PARAMETERS 


TREAD 


0 No initial solution, free stream conditions initialize the flowfield in 
the startup steady-state computations 


1 Initial solution is read from a binary file saved from the previous run 
(default, at the end of each run, solution file is saved as binary) 


2 Initial solution is read as formatted (plot3d form) 
-] Initial solution is read as binary, unsteady motion starts 
-2 Initial solution is read as formatted (plot3d form), unsteady motion 
Starts 


# of timesteps 

Residuals are printed out at every nprint timesteps 

Aerodynamic loads are printed out at every nload timesteps 

Solution variables, q array, are written into "qp.d" file at every delta odvar 


change in unsteady motions (degrees in oscillatory motion, amplitude 
change in plunge) 


Steady state AOA (do not set it to zero, instead set it to 0.0001) 
false N/A 


true _ sinusoidal oscillations in pitch 


false N/A 


_ true straight ramp motion in pitch 


Reduced frequency of the unsteady pitching motion based on the half 
chord, chord length is assumed to be 1. 


Min AOA of the pitching motion 
Max AOA of the pitching motion 
false N/A 

Plunge amplitude in x 

Plunge amplitude in y 


ZF 





PLPHSXD 


PLFREQ 


REYNOLD 


MACH 
VISC 


TURBL 


TRANS 


TIMEACC 


COUR 


NEWTIT 


POTEN 





The Phase angle between x and y amplitudes, in degrees 


The reduced frequency of the plunging motion 
based on the half chord, chord length is assumed to be 1. 


Reynolds number of the freestream flowfield 


Mach number of the freestream flowfield 


false Euler solution 
true Viscous Navier-Stokes solution 


false Laminar flow is assumed. 
true Baldwin-Lomax turbulence model is applied. 
false Fully turbulent flow is assumed. 


true BL transition with Michel criterion is modeled. A steady-state 
solution has to be obtained first. TIMEACC has to be set to true. 


false Variable local time stepping in the computational grid, used in the 
computations of steady state and/or attached flowfields. 


true Constant time stepping everywhere in the computational grid, used 
in the computations of unsteady and separated flowfields 


Courant number of the timestepping (50-1500), determines the timestep of 
the computations based on the minimum grid size and the freestream 
conditions. Its value depends on the computational grid. For diverging 
computations, when the residuals in the output file increases in time, it is 
the sign that its value is to be reduced. 


Number of Newton subiterations in each timestep, applied in unsteady 
flows (2-3), for steady flowfields it is set to 1. 


This variable is reserved for NS-Potential flow interactive solution 
procedure. 


The NS program needs, in general, the following input files: 


ns.in: Input file which defines certain parameters as given above. It is set up for starting 
a generic steady-state solution. For continuing computations or for any changes 
in the input variables you need to edit this file. 


grid.in: A formatted file in which the computational grid coordinates are stored. 


strs.d : A file in which the starting solution is stored. Could be binary or formatted in 
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accordance with the IREAD variable. 


Once you have the input files ready, you can run the NS solver on SGI's by simply 
submitting it with the following command at the prompt: 


run ns 


B. STEADY-STATE CONDITIONS 


C.. 


C.. 


Cc. 


C;. 


Ci 


IREAD, NITER NPRINT, NLOAD ODVAR 

0 3000 10 10 0.01 

ALPHA OSCIL RAMP REDFRE ALFAMND ALFAMXD 
10.0 false false 0.099 0.001 10.0 
PLUNGE PLMX PLMY PLPHSXD PLFREQ 

false 0. 0.10 0. 0.5 

MACH REYNOLD VISC ITURBL TRANS 

0.300 3930000. true ] false 

TIMEACC COUR NEWTIT 


false 100. ] 


C. UNSTEADY CONDITIONS 


Ci. 


C.. 


C.. 


C.. 


C2 


IREAD, NITER NPRINT, NLOAD ODVAR 

-2 13000 20 20 0.01 

ALPHA OSCIL RAMP REDFRE ALFAMND ALFAMXD 
10.0 false false 0.099 0.001 10.0 
PLUNGE PLMX PLMY PLPHSXD PLFREQ 

true 0.0 0.10 0.0 0.5 

MACH REYNOLD VISC ITURBL TRANS 

0.300 3930000. true ] false 

TIMEACC COUR NEWTIT 


true 500. 1 


D. |GRMOVE (MOTION SUBROUTINE) 


subroutine grmove(dalfa,dx,dz) 
include 'coms.f 
dimension xold(nia,nka), zold(nia,nka) 
do i=1,1mx(1) 
do k=1,kmx(1) 
xold(i,k) = xG,k) 
zold(,k) = z(i,k) 


enddo 
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enddo 


if( dalfa .ne. 0.) then 
ca = cos( dalfa ) 
sa =-sin( dalfa ) 
do i=1,imx(1) 
do k=1,kmx(1) 
x(i,k) = xold(i,k) * ca - zold(i,k) * sa 
z(1,k) = zold(i,k) * ca + xold(i,k) * sa 
enddo 
enddo 
endif 
IF( dx .ne. 0. .or. dz .ne. 0. ) then 
if(.not. distort) then 
do k=1,kmx(1) 
do i=1,imx(1) 
x(i,k) = x(i,k) + dx 
z(1,k) = z(i,k) + dz 
enddo 
enddo 
else 
kinner = 45 
kouter = 80 
inner] =89 
yinner2 =283 
iouter] =108 
1outer2 =264 
xrangeinner = 3.5 
xrangeouter = 2. 
yrangeinner = 2. 
yrangeouter = 2. 
do 1=1,imx(1) 
if (i .le. imner1 .and. i .ge. 1) then 
xfactor = 1 


elseif (i .ge. 1inner2 .and. i le. imx(1)) then 


xfactor = 1 


elseif (i lt. iouter] .and. i .gt. iinner1) then 
slope = -(yrangeinner+yrangeouter)/ 


float((iouter1-1)-(inner1+1)) 


xcept = yrangeinner-slope*float(iinner1+1) 
yinput = slope*float(i)+xcept 
xfactor = 0.5*(tanh(yinput)+1.) 


> 


elseif (i .gt. iouter2 .and. 1 It. iinner2) then 
slope = (yrangeinner+yrangeouter)/ 


float((imnner2-1)-(Giouter2+1)) 
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xcept = yrangeinner-slope*float(tinner2+1) 
yinput = slope*float(i+2)+xcept 
xfactor = 0.5*(tanh(yinput)+1.) 


else 
xfactor = 0. 
endif 
do k=1,kmx(1) 
if( k .le. kinner ) then 
factor = 1. 


elseif( k .gt. kinner .and. k .It. kouter ) then 
slope = -(xrangeinner+xrangeouter)/ 
float((kouter-1)-(kinner+1)) 
xcept = xrangeinner-slope*float(kinner+1) 
xinput = slope*float(k)+xcept | 
factor = 0.5*(tanh(xinput)+1.) 
else 
factor = 0. 
endif 
xold(i,k) = x(i,k) 
zold(i,k) = z(i,k). - 
x(i,k) = x(i,k) + dx 
z(1,k) = z(,k) + dz*factor * xfactor 
enddo 
enddo 


endif 
ENDIF 


if( unstdy ) then 


do k=1,kmx(1) 

do 1=1,imx(1) 
xtau(i,k) = (x(i,k) - xold(i,k))/dtau 
ztau(i,k) = (zG,k) - zold(i,k))/dtau 
enddo 


enddo 
else 
do k=1,kmx(1) 
do 1=1,imx(1) 
xtau(i,k) = 0. 
ztau(i,k) = 0. 
enddo 
enddo 
endif 
call metric 


3] 


if( not. unstdy .or. itr .eq. niter ) then 
open(unit=9 1 ,file='grid. out',form="formatted") 
write (91,*) imx(1), kmx(1), twks(1), iwall 

write (91,*) ((x(ik), 1= 1, imx(1)), k = 1,kmx(1) ), 
> ((z(Gi,k), i= 1, mmx(1)), k = 1,kmx(1) ) 
endif 


return 
end 
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APPENDIX B 


A. GRIDGEN INPUT FILE (NACA 4412) 


100 
1.833333 
1.832997 
1.831988 
1.830313 
1.827978 
1.824994 
1.821375 
1.817136 
1.812297 
1.806878 
1.800905 
1.794402 
1.787397 
1.779921 
1.772005 
1.763682 
1.754987 
1.745956 
1.736626 
1.727034 
1.717219 
1.707221 
1.697078 

1.68683 
1.676518 

1.66618 
1.655855 
1.645584 
1.635403 
1.625516 
1.615838 

1.60636 
1.597109 
1.588117 
1.579412 
1.571021 
1.562971 

1.55529 
1.548004 
1.541138 
1.534718 
1.528769 
1.523314 
1.518378 
1.513983 


0 
-1.4E-05 
-5.7E-05 
-0.00013 
-0.00023 
-0.00036 
-0.00051 

-0.0007 
-0.00091 
-0.00116 
-0.00143 
-0.00173 
-0.00206 
-0.00243 
-0.00282 
-0.00325 
-0.00371 
-0.00421 
-0.00473 
-0.00528 
-0.00586 
-0.00646 
-0.00708 
-0.00771 
-0.00834 
-0.00896 
-0.00956 
-0.01013 
-0.01066 
-0.01116 
-0.01167 
-0.01216 
-0.01264 
-0.01307 
-0.01343 
-0.01371 

-0.0139 
-0.01397 
-0.01391 
-0.01371 
-0.01336 
-0.01284 
-0.01216 

-0.0113 
-0.01025 
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1.510148 
1.506894 
1.504238 
1.502195 
1.500779 
1.5 
1.500433 
1.501666 
1.503578 
1.506166 
1.509425 
1.513346 
1.517917 
1.523122 
1.528943 
1.535358 
1.54234 


- 1.549861 


1.557887 
1.566384 
1.575313 
1.584631 
1.594297 
1.604266 
1.614489 
1.624921 

1.63547 
1.645972 
1.656548 
1.667154 
1.677746 
1.688281 
1.698716 
1.709009 

1.71912 
1.729007 
1.738634 
1.747962 
1.756955 

1.76558 
1.773803 
1.781595 
1.788926 
1.795769 
1.802101 
1.807898 
1.813139 
1.817806 
1.821884 
1.825358 
1.828216 


-0.00902 
-0.0076 
-0.00599 
-0.00419 
-0.00219 
0 
0.00471 
0.007159 
0.009655 
0.012185 
0.014728 
0.017264 
0.019767 
0.022211 
0.02457 
0.026815 
0.028917 
0.030851 
0.032591 
0.034114 
0.035397 
0.036424 
0.03718 
0.037653 
0.037837 
0.037726 
0.037324 
0.036695 
0.035879 
0.034885 
0.033726 
0.032416 
0.03097 
0.029403 
0.027731 
0.025973 
0.024145 
0.022265 
0.020352 
0.018425 
0.016503 
0.014604 
0.012748 
0.010954 
0.009241 
0.007626 
0.006128 
0.004762 
0.003544 
0.002489 
0.001607 
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1.830449 0.00091 0 
1.83205 0.000407 0 
1.833012 0.000102 0 
1.833333 07 0 
B. ELLIPTICAL SOLVER PORCESS 
® Select Domain commands from GRIDGEN’s main menu. 





¢ Select the domain to be solved. 
e Set solver attributes (Boundary conditions) 
¢ Foreground control function 
= Enable surface edges (ON) 
> Edge influence (ON) 
= Enter AS delay factor = 10. (Number of gridlines 
from the surface for which these conditions apply.) 
=> Enter angle decay factor = 10. (Angle gridlines 
makes with the boundary itself.) 
¢ Background control function 
= Select Laplace function 
¢ Done with settings 


¢ Run solver. 
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